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Abstract 

Stable, accurate, divergence-free simulation of magnetized supersonic turbulence is a se- 
vere test of numerical MHD schemes and has been surprisingly difficult to achieve due to 
the range of flow conditions present. Here we present a new, higher order-accurate, low 
dissipation numerical method which requires no additional dissipation or local "fixes" 
for stable execution. We describe PPML, a local stencil variant of the popular PPM 
algorithm for solving the equations of compressible ideal magnetohydrodynamics. The 
principal difference between PPML and PPM is that cell interface states are evolved 
rather that reconstructed at every timestep, resulting in a compact stencil. Interface 
states are evolved using Riemann invariants containing all transverse derivative informa- 
tion. The conservation laws are updated in an unsplit fashion, making the scheme fully 
multidimensional. Divergence-free evolution of the magnetic field is maintained using the 
higher order-accurate constrained transport technique of Gardiner and Stone. The accu- 
racy and stability of the scheme is documented against a bank of standard test problems 
drawn from the literature. The method is applied to numerical simulation of supersonic 
MHD turbulence, which is important for many problems in astrophysics, including star 
formation in dark molecular clouds. PPML accurately reproduces in three-dimensions 
a transition to turbulence in highly compressible isothermal gas in a molecular cloud 
model. The low dissipation and wide spectral bandwidth of this method make it an ideal 
candidate for direct turbulence simulations. 

Key words: gas dynamics, magnetohydrodynamics, MHD turbulence, numerical 
methods, PPM, compact stencil 
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1. Introduction 

Piecewise Parabolic Method on a Local Stencil (PPML) 0, Q is a new numerical 
scheme developed for solving multidimensional compressible Euler (HD) and ideal mag- 
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nctohydrodynamic (MHD) equations. The method is based on a piecewise parabolic 
approximation of variables inside individual grid cells. It is third order-accurate in space 
and it implies second order temporal accuracy. This method is an improvement over 
the popular PPM introduced by Colella and Woodward 0, 0] for nonmagnetized flows 
with strong shocks and extended by Dai and Woodward to compressible ideal magne- 
tohydrodynamics PPM has been widely used in computational practice ever since, 
and versions of the PPM gas dynamics scheme have been incorporated into a number of 
codes for astrophysical applications @. 

The main new feature of PPML is the procedure for calculating interface values be- 
tween adjacent cells. Instead of an interpolation procedure used in the original PPM 
formulation, which employs a four-point stencil, PPML relies on the information pre- 
served from a previous time step. The required values are obtained via Riemann invari- 
ants transferred along the characteristic curves of the equations to cell boundaries using 
an approximate parabolic solution within a cell 0, @, @, HI 11 1 ■ To preserve the order 
of the scheme at local extrema, a monotonicity constraint is applied to these interface 
values [Til. [T3I [3, EH • In a multidimensional case a monotonicity preserving method 
from [1 61 ] is additionally applied. The scheme is multidimensional as it keeps terms con- 
taining derivatives with respect to the tangential directions in the equations for wave 
amplitudes. This approach provides the left and right states for the Riemann problem 
based on multidimensional reconstruction. For the ideal MHD case in three dimensions, 
a zero-divergence constraint on the magnetic field is enforced by the use of the Stokes the- 
orem (the so-called constrained transport approach [17|). An updated component of the 
electric field at a cell boundary is calculated by averaging the quantities obtained from 
given components of flux-vectors, taking into account a value of the electric field gradient 
and the information about the sign of the tangential velocity at that boundary [181 ]. 

PPML has been tested on a number HD and MHD problems that demonstrated the 
ability of the method to resolve discontinuous solutions without adding excessive dissi- 
pation. PPML provides a very accurate treatment of strong discontinuities, minimizes 
numerical dissipation of the kinetic and magnetic energy, and substantially improves the 
spectral bandwidth for compressible turbulence models compared to its predecessors [l9| . 

In this paper we present a comprehensive description of the numerical method for 3D 
MHD simulations as well as results of numerical tests. The method has been substantially 
improved and remastered compared to its previous version presented in 0, 0. The 
new PPML features in this paper include: (i) an improved approach to maintaining 
zero divergence for the magnetic field following [ll|, see Section 5; (ii) a monotonicity 



constraint proposed by Rider et al. [14j , as described in Section 6; (iii) an extended 



suite of test problems that includes a comparison with the FLASH3 MHD solver, see 
Section 8; Finally, in Section 9 we illustrate the performance of PPML on a three- 
dimensional problem of highly compressible weakly magnetized forced turbulence with a 
sonic Mach number of 10. The problem of supersonic, super- Alfvenic turbulence with an 
isothermal equation of state proved to be a challenging regime for numerical modeling due 
to the presence of strong rarefactions and very high density contrasts in the flow. PPML 
scheme is perfectly stable numerically on problems of this sort. We also briefly discuss 
the effects of the weak magnetic field on the spectral properties of MHD turbulence and 
obtain good correspondence between our numerical results and observed characteristics 
of supersonic turbulence in star forming molecular clouds. 
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2. PPML description 

Let us consider a homogeneous one-dimensional grid with the spacing h and a function 
q(x,to) = qo{x) defined on this grid at t — to- It is assumed that the function qo{x) can 
be approximated by a parabola inside every grid cell (Fig. [I]): 




Figure 1: Approximation of q(x) inside a cell. 



*>(*) = + f(Ag i + « i (6) (l-0), (i) 

where 

(, = {x-X i _ 1 / 2 )hr 1 , Aqi^q^-q^, 

gf = 6( ?i -l/2(gf + «f)). 
Function satisfies a condition 

q l = hr 1 J q {x)dx. 

Let us consider a linear advection equation 

dF{x,t) 



dt dx 



0, (2) 



where F(x, t) — a q{x, t) and find its solution for the moment t = to + t. On a discrete 
grid we have a number of local Riemann problems which lead to some average states 
q* (x i+ i/2,t) on every interface x i+ i/ 2 between the adjacent cells. The equation ^ has 
only one characteristic defined by dx/dt = a. To find a value <?i+i/2, for example, on the 
right boundary of a cell number i at the time step t = t + r, we suggest to move along 
the characteristic line from the point 2^+1/2 back to the moment t = to and use a value 
from some point on a previous parabola (Fig. [5]) . 
Then for a > we have 

ft+i/afo + r) = q t R (to + r) = 9l L + $ (a© + <zf(l - £)) , (3) 

where 

£ = (a; - a;,-_i/ 2 ) /i -1 = (h- ar) hr 1 = 1 - orAT 1 . (4) 
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All the values in the rhs of (J3]) must be taken from a previous time step t = to . For a < 
a value is defined by a parabola in the cell number i+ 1: 

g i+ l/2(io +T) = q t L +1 (t +t) = q t L +1 + £ [Aq l+1 + gg^l - £)) , 
where £ = -dr/i -1 . 

In monotonic regions, where ft+1/2 € [<& . . . ft+i], it is assumed that q, t R = = 
Qi+1/2 an( i qt — = Qi-i/2- I n non-monotonic regions we must redefine g/' and qt- 
If is a local maximum or minimum, the interpolation function (JTJ) must be constant, 
i.e. qt — qt = qi. If qi is too close to qt or qf", the parabola ([T]) can have an extremum 
inside the grid cell (it happens when |Agj| < |<?.[ 6 ' ) |). In this case we must move this 
extremum to the boundary of the cell. These conditions can be written as 

qt = ft, qt = ft, if (qt - ft) (ft - it) < (5) 

and 

qt = 3g, - 2q«, if Aq t ■ qf > (Aq t ) 2 , 

qt = 3 9l - 2qt , if A Qi ■ qf ] < - (Aq t ) 2 . 
If we know the function q(x), we can compute its average for the interval 
[*i+i/2 - V ■ ■ ■ ^i+1/2] (for y > 0): 

qt+l/2iy)=y~ l j q(x)dx = qt -1/2 yh- 1 

For a > 0, the solution of ([2]) at time i = to + r can be found by averaging over the 
interval [x. l+1/2 - ar . . . x i+1/2 ] , i.e. q*(x i+1 / 2 , to + r) = ft+ 1/2 = ft+i/ 2 ( aT )- For a < 
the zone of influence is [2^+1/2 ■ ■ ■ + a r ] ■ In this case g* (3^+1/2 j £0 + T ) = ft+1/2 = 

vt+1/2 (- flT ), where 

= g£ x + 1/2 y ft" 1 
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Ag 4 -(l-2/3 2/ ft- 1 )g 4 (6) . (7) 



Aq i+1 + (l-2/3yh- 1 )q(% , V > 0. (8) 



The flux on the interface can be computed as 

Fi+1/2 = a + Qi+i/ 2 + a ~<l t +i/2, (9) 

where a + = max(a, 0) = (a + \a\) /2, a~ = min(a, 0) = (a — |a|) /2. We can use an 
arbitrary value for q^ +1 j 2 if a < 0, and for q i 1 ^ 1 ^ 2 if a > 0. 

3. The governing equations 

Let us consider the ideal MHD equations in 3D in the following form: 

d t U + d x F + d y G + d z H = 0. (10) 
Here U is a vector of eight conservative variables, F, G and H arc the fluxes: 

U = (p,pu,pv,pw,B x ,B y ,B z ,E) T , (11) 

F = (pu, pu 2 + p - B^, puv - B x B y , puw - B x B z ,0, 

uB y -vB x ,uB z -wB x ,u(E+p) - B x {uB x + vB y + wB z )) T ', (12) 



G = (pv, puv — B x B y , pv + p — B y , pvw — B y B z , vB x — uB 



0,vB z -wB y ,v(E + p)- B y (uB x +vB y + wB z )) , (13) 

H — (pw, puw — B X B Z , pvw — B y B z , pu? +p — B 2 , wB x — uB z , 

wB y - vB z , 0, w(E +p) - B z (uB x + vB y + wB z )f ', (14) 

where p is the density; u, v and w are the velocity components; B x , B y and B z are the 
magnetic field components; E is the total energy and p is the total pressure: 

BB 

P = P+~- 

We have included the factor l/y/4n in the definition of B. An equation for the total 
energy and an ideal gas equation of state are 

pvv BB 

P= (7-l)pe, 

where 7 - the adiabatic index, e - the specific internal energy. If we denote 
(b x ,b y , b z ) = -1= (B x ,B y , B z ) , b 2 = b 2 x + b 2 y + b 2 z , 



we can write the sound velocity c, Alfven velocity c a , the fast and the slow magneto- 
acoustic velocities c/ iS as 

c = f^- 
V P ' 

C a = \b x \, 



C f,s 



. i( c *+b>)±y(c*+b*) 



b 2 ) -4c 2 b 2 



1/2 



2 v y 2 

We will also deal with a non-conservative form of the MHD equations: 

dtV + AdxV + BdyV + CdyV = 0, 

where 



(15) 
(16) 



V = (p,u,v,w,B x ,B y ,B z ,p) 

The matrices A, B and C can be computed using Jacobians of (|10|) and a transition 
matrix 



M = 



<9U 

Fv' 



For example, the Jacoby matrix A is 
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(17) 



The corresponding eigenvalues are 
\l> 8 = u±c f: X 2 X 7 



U zb C a , A ~' ^ — IliCe, A 



3, 6 



A x ' 8 represent a pair of fast magneto-acoustic waves, A x ' 7 - a pair of Alfven waves, A 
a pair of slow magneto-acoustic waves, X x - an entropy wave, A^ - a magnetic-flux wave. 
The eigenvectors of the Jacobians could have singularities at the points of degeneracy 
of the eigenvalues since the MHD equations are nonstrictly hyperbolic. To avoid those, 
Brio and Wu [2(| suggested a scaled version of the eigenvectors that comes from defining 



(B y ,B z 



By 2 



B' : 



ifB 2 



1 1 



+ Bl ± o, 

otherwise, 



(aj,a s ) = < 



if B 2 + B 2 ^0 or 7 p ^ B 



2 

x ' 



1 1 

—= , —= otherwise. 
V2 V2 



Thus the left and the right eigenvectors are 



' 2Jpc Pyi 2^pc lz ' 2pc< 



rl' s = (paf,±a.f c f ,Ta s c s (3 y sgxiB x ,^a s c s B z sgnB x , 

0,a S T/pcf3y,a sy fpc3 z ,af"/pJ 

1 2,7 _ (n n iko™n §V „„„ u n _l_ §* -^Jk 



. = (0,0,-^^,^-^,0,^,^,0 
■- 7 - ' lUl. sgni?,, ^ sgnS^O, &,=F^f A,, 



1 »■ 6 = ( 0, ± ^ . ± c / Pv s S nB - - ± c f P* s § n ^ . 

^' o /rT^ ' o TtTZ P z ' O „ „ 2 



2^pc y ' 2Jpc '2pc< 



(pa s , ±a s c s ,±a f c f By sgnB Xl ±af c f B z sgnB x , 

0, -ctfy/pcdy, -otf^fpcdz, a s ^pj 



r 3,6 



K = 1 1,0,0,0,0,0,0,--^ 

v\ = (1,0,0, 0,0, 0,0,0) T , 
\\ = (0,0,0,0,1,0,0,0), 
r!j = (0,0,0,0,1,0,0,0) T . 

4. A numerical scheme 

To solve (ITU)) we apply a conservative difference scheme: 

'-'ij.fc ^i.j.fe ^ \ r t+1/2, j, k r i-1/2, j, fc 



Ay V *- 



+ 1/2 _qH+1/2 \ T_ /jjri+1/2 _jjn+l/2 



Ay V *.j+l/2,fe ^"i,j-l/2,kj ^ i,j,fc+l/2 j,j,fc-l/2 



Half-integer indices such as i + 1/2 denote the boundaries of grid cells, half-integer time 
index n + 1/2 means that we use the averaged values of the fluxes over a time step t in 
order to get a second-order temporal accuracy. 

The solution inside every grid cell is approximated by a parabola along any Cartesian 
grid axis. The boundary values for each parabola are determined from a conservation 
property of Riemann invariants that remain constant along the characteristics of the ini- 
tial linearized system of equations. Parabolae must be built using the primitive variables 
P^|) . so we need to consider a non-conservative form of the MHD equations (TTSJ). 

We can expand A, B and C in (fT5|) into their eigenvectors. For example, for the 
x-direction: 

A = R x K x L Xl (19) 

where R x is a matrix with columns filled by the right eigenvectors (p = 1, ... ,8), 
L x = R^ 1 is an inverse matrix, with rows filled by the left eigenvectors l p . A x is a 
diagonal matrix of the eigenvalues: (A. x )ij — for i ^ j, (A x )jj — X p for i = j = p. 

To construct piecewise parabolae for every time step, one needs to define the states 
on the cell edges and the states at their centers. For simplicity, let us further consider a 
ID case: 

d t V(x,t)+Ad x V(x,t)=0. (20) 
Inserting (|19[) into (j2"U)) and multiplying by the L matrix from the left, we arrive at 

Ld t V + ALd x V = 0. (21) 

Let us expand a vector V(x, i) into the local basis of the right eigenvectors r p , which are 
fixed in every cell: 

V(:M) = 5> p (a:,i)i- p . (22) 
p 



Inserting (|22[) into (121|) , we arrive at 

d t a p + X p d x a p = 0, p = l,...,8. (23) 

The equations (|23|) mean that the coefficients a p (x,t) in the expansion (f2"2"|) (the wave 
amplitudes) must be constant along the characteristics x p (t): 

dx p _ 

i.e. a p (x, t) are Riemann invariants. A value of Riemann invariant on the boundary of a 
cell (x = x i+ i/2) at the moment t + t could be computed using its value at the moment 
t as 

a p (x i+1/2 ,t + r)=a p (x l+1/2 -X p r, t). (24) 



Fig. GUshows two adjacent cells i and i + 1. The characteristics in the cell i have index 
Pi, in the cell i + 1 - index p2- One of the characteristics x Pl (£) with the eigenvalue 
A Pl > is shown in the cell i, another one x P2 (t) with the eigenvalue X P2 < is shown 
in the cell i + According to the amplitude of a wave at point 3, which propagates 
inside the cell i along the characteristic x Pl (t) with the eigenvalue A Pl , is equal to its 
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value at point 1. In the same way the amplitude of a wave at point 3, which propagates 
inside the cell i + 1 along the characteristic x P2 (t), is equal to its value at point 2. 

The state at point 3, which is computed according to (|22|) by summation with respect 
to all the eigenvectors, fixed in cell i, with A Pl > will be on the left side of the interface. 
Let V L denote this value. Similarly, let denote the state at point 3 on the right side 
of the interface which is computed by summation with respect to all the eigenvectors 
fixed in cell i + 1, with A P2 < 0. 

The amplitudes of waves a p (x p ,t) at the point x p — 3;,-+ 1/2 — A p r in jM]), which 
influence the right boundary of cell i (A p > 0), can be computed by multiplying the 
expansion (|22p by the left eigenvectors, fixed in cell i: 

a p (x p ,t) = l p V(x p ,t), A p >0, (25) 

where 

V(x p ) ee / V(x)dx. (26) 

Xi+1/2 - x p J x p 

We can use arbitrary values for the wave amplitudes for A p < because these waves 
have no effect on the right boundary of the cell and will be omitted in the sum (|22|) . For 
convenience let them be 

a p (x p ,t) = l p V(xi,t), \ p <0. (27) 
With the help of 0-function we can rearrange (|23)) - ([2T|) as 

a p (x p , t) = e(X p ) (l p V(x p , t)) + (1 - 9(A p )) (l* V(ii, t)) . (28) 

Multiplying (J2HJ) by r p and summing over all p such as A p > according to (f2"2")l , after 
some simple manipulations we obtain the boundary value at time t + r: 



V L (x t+1/2 ,t + r)=V(x l ,t)+ r p [l p [V(x p ,t)-V(x l ,t)j\. (29) 

p (Ap>0) 

If we consider cell i + 1 and waves with A p < wc will obtain a similar expression 
for the value V fl on the right side of the interface at time t + r: 

V R (x l+1/2 ,t + T)=V(x l+1 ,t)+ J2 rP [l p (v(z p ,i)- Vfc+i.t))] . (30) 

p (Af<0) 
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The left and the right eigenvectors in (f29 |) - ([30|) are fixed in every cell. To compute 
them we can use a state from any point x inside the cell - it has been shown by numerical 
experiments that this choice has no influence on the solution. We suggest using the values 
of states in the centres of the cells, i.e. l p = l p (V(xi, t)). 




Figure 4: Approximation of V(ai) in the adjacent cells. 



A possible approximation of the component V(x) of the vector-function V(ir, t) inside 
grid cells i and i + 1 at some time step is shown in Fig. 01 The arrows point to the values 
on the left and the right sides of the interface, note that V L ^ V R . The dotted lines are 
the average values of V(x): 



Vt = 



V(x)dx. 



-1/2 



In order to solve the 3D problem we split the initial set (fT5)) by the space variables and 
solve the ID equations separately for the x-, y- and z-directions. However, in this case 
we have an additional change in the quantities because of the fluxes in the orthogonal 
directions. For example, a flux in the y-direction will affect the quantities at point 1 
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Figure 5: A 2D mesh. 



between cells and (i + (see Fig. To obtain the correct result, we can solve 
equation (Tl3| in the x-direction considering the terms B d y Y and C d z ~V as sources. 
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Then instead of (f23|) we arrive at 



d t a p + \ p d x a p = -D p , p = l,...,8, (31) 
where D p are the components of the vector 

T> = L x (Bd y V)+L x {Cd z V). (32) 
The components of the partial derivatives in (f3"2"| can be calculated as 

V p - y p 

^ yp _ i,j + l/2,k ^ i,j-l/2,k ^gg^ 

The solution of (f3"Tj) can be obtained from a Taylor series expansion of a p (x,t) near the 
boundary x = x i+ i/ 2 . It is similar to (|24p but has the additional term 

a p {x t+l/2 ,t + T) =a p {x t+1/2 -\ p T,t)-D p {t) T -. (34) 

Then instead of (f2"5)) - (|30[) we arrive at 



V L (x i+1/2 ,t + T)=V{x i ,t) + r p [l p {y(x p ,t)-V(x h t)-D p (t)^\, (35) 

p (Ap>0) 

V R (x l+1/2 ,t + r)=V(x l+1 ,t)+ r p [l p (v(x p ,t)^V(x l+1 ,t)~D p (t)^j' . 

p (Ap<0) 

(36) 

In f3"4"|) - (|3"6"|) we have omitted indices j and k. The state at point 2 on Fig. [5] is computed 
in a similar way considering the x-derivative as a source. 

To obtain states V. i+1 / 2 in PPML, we solve the Riemann problem for the V L and 
V R states on every interface using, e.g., the Roe solver [2l[ or the HLLD solver [23]: 

V 4+1/2 = i?(V L ,V*), (37) 

where R is the Riemann solver. Note, that in the original PPM, the states V i+1 / 2 
are obtained through monotonic interpolation We then apply a monotonicity- and 
extrema-preserving procedure proposed by Rider et al. [3] to the values of Vj +1 / 2 , as 
described in Section[SJ Finally, we modify the resulting interface states with the standard 
PPM monotonicity procedure ©-(d]). 

So far we have obtained the boundary values of piecewise parabolae at time t + r in 
each grid cell. Now we need to define the fluxes to compute the new central states. Again 
we can use the Roe solver [HI or the HLLD solver j22| with V L and ~V R from f35 ]) -([36 ]l 
but in this case the numerical scheme will have the first order of temporal accuracy. To 
design a second order scheme, we must average the amplitudes a p (x,t) over the zones 
of influence. 

Fig. El shows a set of characteristics corresponding to waves with \ p > 0. The 
characteristic x 1 (t) has the maximum eigenvalue A 1 . Point 1 is the point of intersection 
between this characteristic and the piecewise parabola at time t. Obviously, only the 
zone between the interface x = x i+ i/ 2 and point 1 affects the left boundary state at 
point 2. 
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Figure 6: Characteristics inside a grid cell. 



If we consider a wave which propagates along the characteristic with X p > inside 
cell i, its averaged amplitude on the interface x = x i+1 / 2 at time t + t can be calculated 
as 

"f+i/2 = ^ / <* P (x)dx, A*>0. (38) 

Multiplying the expansion (|22|) by the left eigenvectors, fixed in cell i, will yield 

a p (x) =P V(x), A p >0. (39) 
Inserting (|39[) into (|3"5)l and removing the factor \ p from the integral, we arrive at 

— p L.p 

a i+l/2 = 1 Vi+1/2 > 



V i+l/2 



After that we can arrive at 
— l — L,l 

V =v. i+1/2 



\ p T 



V(x)dx, X p >0. 



a=i+l/2-A p r 



E : 

p (Ap>0) 



— L. p — L, 1 <t- 

V , — V , — D p — 

V i+l/2 V 1+1/2 u 2 



(40) 



(41) 



which is similar to ([35)1 . Here V i+1 / 2 is the averaged by formula (|40|) solution V(x) at 
time £ over the zone of influence of the wave in cell i, corresponding to the maximum 
eigenvalue (A 1 > 0). 

For cell i + 1 and the negative eigenvalues we arrive at 



— R — R, l 

v = v i+1/2 



„ r /—r,p —r, 
E rP 1P ( v m/2-V l+1 



p (Ap<0) 



1 — R,p — R, l r 
/2 2 



(42) 



V <+1 



+1 / 2 ~ |Af|r 



V(s)dx, A p <0. 



(43) 
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i.i+ 



zi,j+1/2 v 



■Zi+1,j+1/2 



^Zi+1/2,j+1/2 

E z i+1/2j 



Figure 7: Computation of the component E z at a node. 



To get the second order of accuracy for time we must compute fluxes on every interface 

— l — R 

solving the Riemann problem with V and V from (|4i p - (|42|) . Note that the values of 

components of the integrals (jit!)) and (j4"3"|) can be computed as ([7]) and respectively. 

The time step t in (JTSJ) is obtained from the Courant condition 

. { Ax Ay Az ) fJJ . 

T = (Tmin < : : , : : t. , : : > . (44) 

HJ.fc [ \Ui,j,k\ + Cf ijjJe \V id>k \ + C y flyj k \WiJ, k \ + C] l j k J 

where a is the Courant number and Cy, A, and c| are the fast magneto- acoustic velocities 
along the coordinate directions. 

5. Zero divergence constraint for the magnetic field 

Our numerical method must provide numerical solutions that satisfy a condition 

divB = 0. (45) 



There are several approaches to this problem in the literature [181 1231 . [24j, |25|, [17|, [26| . 
In our numerical scheme we used an unsplit Godunov method for ideal MHD with a 
constrained transport developed in that is based on the Stokes theorem 

^=-VxE. (46) 

Calculated from (|46p the magnetic field obviously satisfies condition (|43|) . To apply the 
Stokes theorem in 3D case, we define the magnetic field components on the cell faces 
and the electric field components on the cell edges. The algorithm exploits the fact that 
some components of the fluxes F, G and H are actually the components of the electric 
field. 
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For simplicity let us consider a 2D case. Thus the sixth component of the flux F is the 
z-component of the electric field reversed in sign, which correspond to the left and right 
boundaries of a cell. The fifth component of the flux G is the z-component of electric 
field, which correspond to the top and bottom boundaries of a cell. Computing Taylor 
series for these components near the nodes of the mesh and averaging them, we obtain 
the components of the electric field at the nodes (Fig. [7]) . 



771 _ ( j^left . T^right 

&z J+1/2J + 1/2 - ^ J+ l/2j + l/2 + J+1/2J + 1/2 



. rjitop . rpbottom \ 

+D z 1+1/2 j + 1/2 + i+l/2j + l/2 J > 



where for example 



E 



left 



i+l/2j + l/2 - E * iJ+1/2 "I" q x 

fdE : 



dE z 



dE z 



Ox 



JJ + l/2 



\ dx 
dE. 



ij+1/2 



dE z 
dx 



( dE z 
V dx 



Ax 

v i,j + l/2 > 0, 
v i,j + l/2 < 0, 

1 v i,j + l/2 = 0. 



These values of E z i+i/2j+i/2 are use d in a discrete version of (^46]) : 



B 



'Tl-t-l r>n (771 \ 

x i+l/2, 3 - °x i+1/2 ,j - i+l/2J+l/2 - & z i+ l/2j-l/2 ) , 

B ytj + l/2 = B y i,j + l/2 + ^ ( E z i+l/2,j+l/2 - E z i-l/2,j+l/2j ■ 

Updated components of magnetic field at the center of (z,j)-cell are computed by 
averaging: 



b t :+] = l( B n+l 



B 



n+l 



j 2 V x '-VSJ a i+l/2,j) ' 

i,j 2 V v 1 ~ r a i,i+i/2 ,1 ' 

The magnetic field computed this way will automatically satisfy (|45[1 . To make sure 
that this is indeed the case, divB should be approximated by the following expression 



(divB) i+1//2j+1 / 2 — YXx \ x + Bx - Bx hj ~ Bx 
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6. Monotonicity constraints 



A standard PPM monotonicity preserving procedure ([S])-© is insufficient for the 
ideal MHD case. We need to use additional procedures to suppress spurious oscillations. 

Procedure 1. To keep a solution monotonic without reducing the order of the scheme 
and to preserve all the local extrema in MHD simulations we can follow a number of ap- 
proaches described in the literature 12, HI 14 , HI ■ In our method we rely on a Piecewise 



Parabolic Accurate Monotonicity- and Extrema-Preserving procedure described in 14 1. 

We compute the wave amplitudes a? for the central states V$ and a^ ±1 ^ 2 for the 
interface states ~Vi±\/2- We then calculate new values for wave amplitudes of the interface 
states: 



( i ■' — 1 1 1 1 ' U ] 1 -.1 1 i I i \ ' I \ ' I 1 ■' 

and 



a' i±1/2 = median ( of , aj ±1/2 , ar j±1 



a 



= median ^af, CKy_i 3<i'' - >n ' 



i±l/2 mcuJOU ^"i ' "iil/2' ^ u i^l/2 

where the median function is defined in a usual way 

median (a, b, c) = a + minmod (6 — a, c — a) 
through the minmod function 

minmod (a, fe) = — (sign(a) + sign(&)) min (|a|, |&|) . 

If 

a f±i/2 — a i±i/2 au Pi P roce dure is completed. Otherwise, we compute a set of 
fifth-order WENO interface values o^ ±] y 2 # , using Algorithm 2.2.4 from [3], and check 
for a local extremum. If = a i ^ or au we a PPly Algorithm 2.1.2, 4(b), otherwise 

the region is monotonic but too steep to be approximated with the values obtained from 
the Riemann solver ([57]) and we apply Algorithm 2.1.2, 4(c). 

Procedure 2. To keep a solution monotonic in a multidimensional case, we employ 
a method proposed in [16j ]. 

Let us consider a 2D case for simplicity. A parabola that approximates a solution 
along the ir-axis for every component V(x) of a state ~V(x,t) at some point in time can 
be defined as 



V(x) = V t j + <t>{V) 

where 



(x - Xl ) + ^i((x- Xl f - 



__ j4+l/2J ~ ^-l/2j _ _ Vj+l/2.j — 2Vj,j + ^-l/2,j 

Sl - J ~ ~Kx ' ai ' j 

As a limiting function <p(V) we can use that described in flit ]: 

4>(V) = min f 1 |^-max(^ m )| 



Vij - max(Vi_i/ 2i j, V i+1 / 2> j, V it j_i/ 2 , V ij+1 / 2 )\ ' 
\Vi,i - min(Vj jTra )| 



- min(V^_ 1 / 2 j, Vi +1 / 2j , I^j_i/2, ^j+1/2) 



(47) 



where I = i — 2, i — 1, i, i + 1, i + 2, m = j — 2, j — + 1, j + 2 except (Z, m) = (i, j). 
In a 3D case, the limiting function (|4"7f must include all the neighbors of the cell 
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7. A FORTRAN implementation 



The algorithm for computing the left and the right boundary values V L and 
35j)d36j) could be implemented in FORTRAN this way: 



integer n, n2 
real dt , dx 

real VL(8), VC(8), VR(8) , Vm(8) , V(8), D(8), dVy(8) 
real Lambda (8) , sumL(8), VLnew(8) , VLinterf ace (8) 
real B(8,8), L(8,8), R(8,8) 



call Eigenvalues (VC, Lambda) 
if (Lambda(l) .gt.O.) then 



VLnew=0 . 

xi=l . - Lambda(l)*dt/dx 
call Vxi(xi,VL,VC,VR,Vm) 



compute the eigenvalues at the center of the 
cell i (VC - a state in the center) 
if A 1 > then compute a new state VLnew 
on the left side of the interface between the 
cells i and i + 1 

£ in CO) for A 1 (see Q) 
formula (J); VL, VR - the left and the right 
boundary values of a parabola, Vm - the result 
compute the matrix B 

compute the left (L) and right (R) eigenvectors 



call MatrixB(VC,B) 
call Vectors(VC,L,R) 
D=0. 

do n=l,8 
do n2=l,8 

D(n)=D(n)+B(n,n2)*dVy(n2)/dy ! (BP d y V p ) (see fl5J)-(|55])) 
enddo 
enddo 
D=D*dt/2. 
sumL=0 . 
do n=l,8 

if (Lambda (n) . gt . . ) then ! only these waves affect the interface 
xi=l.- Lambda (n)*dt/dx ! £ in Q for A p 

call Vxi(xi,VL,VC,VR,V) ! V - the result (a state at the point £) 
do n2=l,8 

sumL(n)=sumL(n)+L(n,n2)*(V(n2)-Vm(n2)-D(n2)) ! apart of ([55 
enddo 
do n2=l,8 

VLnew(n2)=VLnew(n2)+R(n2,n)*sumL(n) ! a part of ([55]) 
enddo 



end if 
enddo 

VLnew=Vm+VLnew 
else 

VLnew= VLinterface 
end if 



! the result for (1551) 

! if the maximum eigenvalue A 1 < 

! we keep the old value on the left side of the interface 
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A FORTRAN code for computing V and V is similar. We only need to replace 
the function call call Vxi (xi , . . . ) with one that computes the integral |Q or JHJ), using 
xi=Lambda ( 1 ) *dt/dx 
or 

xi=-Lambda(8) *dt/dx , 

respectively. Lambda (8) is the maximum absolute value of the negative eigenvalue. 

To obtain a monotonic solution for a 3D MHD problem, we suggest using the following 
algorithm: 

— L — R 

1. Compute the average interface states V and V from (|4Tjl - lj42)l . 

— l — R 

2. Solve Riemann problem between V and V to determine the fluxes. 

3. Use a conservative difference scheme (|18p to compute new central states. 

4. Modify the magnetic field at the centers accordingly, see section [5] 

5. Compute the interface states V L and ~V R from (f3"5|) - ([3"6")) . 

6. Solve Riemann problem between ~V L and V R to determine the new interface states. 

7. Apply procedure 1 in ^-direction (section E]). 

8. Apply procedure 1 in y-direction (section E]). 

9. Apply procedure 1 in z-direction (section [6|) . 

10. Apply procedure 2 (section [5]). 

11. Apply PPM procedure ©-©■ 

Note that only the interface states are modified by the monotonicity preserving pro- 
cedures. The solution of the Riemann problem joins the interface values, but the mono- 
tonicity procedures split them again into the left and right values. 



8. Numerical tests 

8.1. Riemann problem with multiple weak discontinuities 



This is a ID problem from [27J. An interval x G [0 ... 1] is divided in two by x = 0.5. 



The left and the right states at the initial moment are defined as 

(p L ,u L ,v L ,w L ,B^,B z L ,p L ) = (1.08,1.2,0.01,0.5,3.6,2,0.95), 

(p R , u R , v R , w R , B R , B R ,p R ) = (1, 0, 0, 0, 4, 2, 1), 

B x = 2, 7 = 5/3, N = 512. The solution involves two fast shocks with Mach numbers 
1.22 and 1.28, two slow shocks with Mach numbers 1.09 and 1.07, two rotational and 
one contact discontinuities. The solution for the moment t = 0.2 is presented in Figs. [8} 
1101 The solid line represents the exact solution and points represent the numerical one. 
PPML produces very sharp fronts resolved with only a few grid points. 
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Figure 8: Riemann problem with multiple weak discontinuities. Density and pressure distributions. 



1 - 

o 



Figure 9: Same as in Fig. |8]but for y- and ^-components of the magnetic field. 

V 



Figure 10: Same as in Fig. [8] but for x- and y-components of the velocity. 
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8.2. Numerical dissipation and decay of Alfven waves 

Numerical calculations on a discrete grid always lead to energy loss due to numerical 
dissipation. In order to estimate the properties of numerical dissipation of the PPML 
ideal MHD scheme, we used the test problem from [28{ and followed a decay of two- 
dimensional Alfven wave. We used a standing wave propagating along the grid diagonal 
with initial conditions 

Sv x = v amp c a sin(k x x + k y y), 

5p = Sp = 8v x = 5v y = SB X = SB y — SB Z = 

in a stationary background flow with po = 1, po = 1, B x — 1, B y = B z = 0. This gives 
the sound speed c = 1.291 and the Alfven velocity c a — 0.7071. The computational 
domain is a square box with size L = 1 divided into 64 x 64 grid cells. The wavenumbers 

k x = k y — 2n/L, the total wave number k — ^jk\ + k^ = \/2(2w/L), the initial peak 

amplitude v amp = 0.1, adiabatic exponent 7 = 5/3. Computations were carried out with 
a Courant number a — 0.4. We used the periodic boundary conditions. 

Figure [Tl] shows the envelope for the maxima of z-component of the magnetic field 
and velocity obtained with PPML and PPM reconstruction procedures as a function of 
time. While both schemes show very low dissipation, PPML dissipation is even smaller 
than that of PPM. 

8.3. Travelling circularly polarized Alfven wave problem 

This problem was suggested in [23| as a test for numerical accuracy of smooth flow 
solutions. The circularly polarized Alfven wave propagates at an angle of a — 30° with 
respect to an axis x in the domain [0, 1/cosa] x [0, 1/sina]. The initial conditions are 

p=l, U|| = 0, vj_ = 0.1 sin(2<), w = 0.1 cos(2?r£) 

S||=l, B ± = 0.1sin(27r£), B z = 0.1 cos(27r£), p = 0.1, 

where £ = a; cos a + ysina. For convenience the parallel and the orthogonal to the 
direction of Alfven wave propagation components of the velocity and the magnetic field 
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Figure 12: The orthogonal component m travelling Alfven wave in computations on meshes with 
N = 8, 16, 32 and 64. 

are used instead the of the components u, v, B x B y . For example B\\ = S^cosa + 
By sin a, B± = B y cos a — B x sin a. Alfven wave travels to the point (x, y) = (0, 0) with 
the velocity B\\/yfp = 1. Note that the wave becomes standing if vu = 1. 

The problem was solved on a set of rectangular N x 2N meshes with N — 8, 16, 32 
and 64. The averaged relative numerical errors were estimated as 

N 2N 

W) = * X V 2JV . ioi-U = v ± ,w,B ± ,B z , (48) 

E E \u%\ 

where the solution on the mesh iV = 128 regarded as the exact one Ufy The rate of 
convergence was calculated as follows 

R N = log 2 (S N/2 /S N ) , (49) 



Table 1: Travelling Alfven wave. The average relative errors and the rates of convergence at t = 0.5. 









N 


S N 


Rn 


8 


2.2384 x 1CT 1 




16 


5.7258 x 10" 2 


1.967 


32 


1.7031 x 10~ 2 


1.755 


64 


5.0365 x 10~ 3 


1.771 
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Figure 13: Rotor problem. The contours represent thirty levels of the density in the range from 1.3 to 
13.5 and the pressure in the range from 0.12 to 2.1 at t = 0.15. 




Figure 14: Rotor problem. The lines represent thirty levels of absolute value of the magnetic field in the 
range from 0.32 to 2.288 and the Mach numbers in the range between 0.144 and 4.27 at t = 0.15. 



where S n is an averaged value: 

Sn = \ (S N {v ± ) + S N (w) + 5 N (B ± ) + 5 N (B Z )) . (50) 

The calculations were carried out up to t = 5 with a Courant number a = 0.4 and 7 = 
5/3. We applied periodic boundary conditions. Figure [T2l demonstrates the convergence 
of the numerical solution. It shows the orthogonal component B± in travelling Alfven 
wave for computations on meshes of different size N. Table [T] gives average relative 
numerical errors and rates of convergence obtained for the PPML scheme. 
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8.4- Rotor problem 

The Rotor problem was suggested in [26| and has been widely used to test numerical 
schemes [e.g., 23, IH, 3(|. It turned out to be a hard nut to crack for many codes due to 
the appearance of negative pressure values [iij ]. 

The computational domain in this case is a square [0, 1] x [0, 1] with a uniform pressure 
p = 1 and magnetic field components B x — 5/-\/47r, B y = 0. There is a rotating disk 
of dense fluid at the center with a radius ro = 0.1. For r < ro we specify p = 10, 
u = —vo(y — 0.5)/ro, v = vq(x — 0.5)/ro, where r — y/ (x — 0.5) 2 + (y — 0.5) 2 , vq = 2. 
For r > r% = 0.115 the fluid is initially at rest (u = v = 0) with density p = 1. In the 
intermediate zone tq < r < r\ we use a linear interpolation of the variables: p = 1 + 9/, 
it = —fvo(y — 0.5)/r, u = Jvq(x — 0.5)/r, f = (r% — r)/(ri — ro). In this setup, the 
initial configuration is unbalanced due to centrifugal forces. The rotating fluid will tend 
to equilibrate, while the magnetic field holds the oblate shape of the rotor. 

The computations were carried out on a set of N x N meshes with N = 50, 100, 200 
and 400, with a Courant number a — 0.4 and 7 = 1.4 until t — 0.15. The boundary 
conditions are obtained through zero-order interpolation. Figures Q1J] and [TJ] show the 
flow fields for N = 400. 



Table 2: The rotor problem. The average relative errors and the rates of convergence in the numerical 
codes at t = 0.15. 



N 


PPML 


Flash3 USM-MEC 


<5jv 


Rn 


S N 


Rn 


50 


9.4274 x 10~ 2 




1.1470 x 10" 1 




100 


4.5204 x 10" 2 


1.06 


5.9800 x 10~ 2 


0.94 


200 


1.9262 x 10~ 2 


1.24 


2.5000 x 10~ 2 


1.26 



We compare numerical solutions obtained with PPML with those presented in 23, 2^, 
[30j . In Table[5]the average relative numerical errors and the rates of convergence are given 
for PPML and for Flash3 USM-MEC (unsplit staggered mesh algorithm with modified 
electric field construction introduced in Both PPML and Flash3 USM-MEC codes 

use a Roe solver [2l| for this test. The relative numerical errors were computed using 
eq. (j4"8|) , where for Ufj we used the highest resolution result (N = 400). The average 
error 5n(U) is defined as the average Sn(U) for all non-zero variables U. The rate 
of convergence is estimated as in (|4T))) . PPML results are more accurate and have a 
comparable rate of (self-)convergence with those from the new Flash3 MHD solver. 

8.5. Orszag-Tang vortex problem 



This problem was suggested in 31[ and since then has been used in many papers as 
a standard test problem for numerical codes in 2D MHD. It involves formation and an 
interaction of multiple shocks and a transition to supersonic turbulence. 

In the computational domain [0, 1] x [0, 1], we set a uniform density p — 25/(367r) and 
pressure p = 5/(127r) with 7 = 5/3 (in this case the sound velocity c = \fjpjp =1). The 
initial velocities and components of the magnetic field are set using harmonic functions: 
11 = — sin27ry, v = sin27ra;, w = 0, B x = — i?osin27ry, B y — Bosm4-irx, B z — 0, where 
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Bq = 1/v4tt. Despite such smooth initial conditions the fluid motion becomes very 
complex. 

We carried out computations on N x N meshes with several values of N using periodic 
boundary conditions and a Courant number a = 0.3. Figure fTBl demonstrates the pressure 
distribution at time t = 0.5 for N = 256. In Fig. [TC] the pressure distributions along the 
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lines y = 0.3125 (j — 83) and y — 0.4277 (j = 112) are shown to illustrate the accuracy 
and sharpness of the main flow features. 

Table [3] contains the average relative numerical errors and the rates of convergence 
for PPML and Flash3 USM-MEC (3^ solvers at t = 0.5. As an "exact" solution U^, we 
used a solution obtained on the grid with N = 400. The corresponding PPML results 
are more accurate and demonstrate better convergence. 

Table 3: Orszag-Tang vortex problem. The averaged relative errors and the rates of convergence in the 
numerical codes at t = 0.5. 



N 


PPML 


Flash3 USM-MEC 


6 N 


i?jV 


8 N 


Rn 


50 


8.9095 x 10~ 2 




1.0160 x 10" 1 




100 


4.4249 x 10~ 2 


1.013 


5.2200 x 10" 2 


0.964 


200 


1.8851 x 10~ 2 


1.235 


1.9900 x 10~ 2 


1.390 



9. A compressible turbulence simulation 

Magnetized supersonic turbulence plays an important role in statistical star formation 



theories [35j. This stimulated development of accurate numerical methods suitable for 
modeling turbulent molecular clouds. One of the motivations behind the design of PPML 
has been a need for an MHD scheme with low numerical dissipation comparable or better 
than that of PPM. In this section we illustrate the performance of PPML on a challenging 
problem of forced super- Alfvenic turbulence. Some numerical methods that successfully 
pass the tests discussed above turn unstable on this application. Since adding more 
dissipation where needed - the usual way to cure for "blow ups" caused by numerical 
instabilities - would ultimately damage the derived statistics of turbulence [33| , the issue 
of inherent stability of numerical methods is crucial for both supersonic turbulence and 
star formation simulations. 

For illustrative purposes, we present here a simulation of weakly magnetized super- 
sonic turbulence. In this experiment, turbulence in a periodic domain of linear dimension 
L = 1 is driven by a large-scale solenoidal force for 8 flow-crossing times ta = L/2M S . 
At time t = 0, a uniform gas with density p = 1 is permeated by a weak uniform mag- 
netic field Bo || x, such that (3q = 2p/Bq = 20. We apply an initial large-scale velocity 
field that corresponds to an rms sonic Mach number M s ~ 10 and assume an isothermal 
equation of state (c = 1) to mimic the average physical conditions in the dense parts of 
molecular clouds (n = 10 3 cm -3 , T = 10 K). 

The evolution begins with the formation of strong shocks on "caustics" of the initial 
velocity field. Shock interactions then cascade the initial kinetic energy of large-scale 
motion of the gas down to smaller and smaller scales a la Kolmogorov-Richardson, see 
Fig. [T71 The magnetic field gets amplified by a factor of about 50 via the small-scale 
dynamo action [34[. The large-scale solenoidal force (acceleration) keeps the rms sonic 
Mach number roughly constant at M s w 10. The evolution of kinetic and magnetic 
energies is shown in Fig. [THJ left panel. Also shown is max(|V-B|) as a function of 
time during this simulation. The method keeps the absolute value of the divergence 
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Figure 17: Supersonic turbulence simulation with PPML on a 512 3 grid. Four snapshots show the 
density field on a slice x = illustrating a transition to fully developed turbulence with M s as Ma ~ 10. 
The transition includes formation of first strong shocks on the caustics of the initial solenoidal velocity 
field (t = 0.5^, top-left), shock interactions and active development of first shear instabilities (t = lt^, 
top-right), and a gradual transition to a statistical steady state (t = 2t^ and At^, bottom row). The 
standard logarithmic black-red-white color ramp shows high-density regions in light-red and rarefactions 
in black. 

of magnetic field below 10~ 12 at all times, even after 70,000 integration time steps (if 
double precision is used). After about 4 crossing times of evolution, the system completes 
a transition to a fully developed isotropic state with M s « 10 and Ma « 10. The right 
panel of Fig. [T8l illustrates spectral characteristics of this saturated state by showing the 
time-average (over 25 flow snapshots taken between t — 4t<j and t = 8td) power spectra 
for the density, velocity, and magnetic field strength. 

The velocity spectrum has an extended scaling range with a slope of about —2, as in 
the Burgers turbulence, similar to the corresponding scaling in non-magnetized flows at 
high Mach numbers 32]. This is expected, as turbulence here is only weakly magnetized. 
The density spectrum slope of about —0.7 is again consistent with our previous results 
for non-magnetized flows obtained with PPM (—1.0 at M s = 6) and with an anticipated 
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Figure 18: Supersonic turbulence simulation with PPML on a 512 3 grid. Time-evolution of kinetic and 
magnetic energy and maximum absolute value of V-B (left panel) and turbulent power spectra for the 
velocity, density and magnetic field (right panel). 



trend towards a flat "white noise" spectrum at M s = oo. The magnetic energy spectrum 
does not show a clear scaling range, as expected at this modest resolution, assuming the 
effective magnetic Prandtl number of PPML is of order unity. We also looked at more 
advanced spectral characteristics for compressible flows, such as the power spectrum of 
p 1 / 3 v, and found a slope of —1.7. This power spectrum related to the energy transfer 
rate in wavenumber space is insensitive to the turbulent Mach number and should have 
a Kolmogorov —5/3 slope [32| in both incompressible and highly compressible regimes, 
although a steeper scaling does occur due to intermittency [36J. This is true for both 
non-magnetized and weakly magnetized flows. 

We have also carried out two additional simulations of the same kind but with higher 
degrees of magnetization, (3q — 2 and 0.2 37|. While the saturated turbulent state 
in the (3q = 2 variant is still super- Alfvenic with Ma ~ 3 and the magnetic energy is 
about 3 times smaller than the kinetic energy, the trans- Alfvenic case, (3q — 0.2, reaches 
an equipartition of kinetic and magnetic energies. In both cases, PPML proved to be 
perfectly stable at a Courant number a = 0.2, as in the super- Alfvenic case (3q = 20 
discussed above. 

Our approach to handle the stability issues in MHD turbulence simulations with 
PPML is as follows: (i) we use locally multidimensional reconstruction that improves 
the quality of Right and Left interface states and helps to avoid numerous well-known 
pathologies, such as "carbuncles", etc. (38j; (ii) we control the quality of these states 
before moving forward with the flux calculation; (iii) if we find that the states are not 
satisfactory, we reduce the order of reconstruction to linear or further skip the whole 
reconstruction step; (iv) we use only nonlinear Riemann solver in all cases lacking the 
reconstruction step. 

Overall, the derived spectral properties of weakly magnetized highly compressible 
turbulence demonstrate that low dissipation and wide spectral bandwidth of PPML 
make it an ideal numerical scheme for large-scale simulations of magnetized supersonic 
turbulence. 
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10. Conclusions 

In this paper we presented PPML, a new numerical method for compressible ideal 
MHD that is based on the piecewise parabolic approximation. Interface values for the 
interpolation parabolae in every grid cell are defined with the help of the Riemann in- 
variants which remain constant along the characteristics. The monotonicity of the states 
on the interfaces between adjacent cells is provided by the monotonicity- and extrema- 
preserving procedure from [14j |. The scheme is fully multidimensional as it includes the 
terms corresponding to the tangential directions in the amplitude equations. This helps 
to avoid numerous well-known pathologies, such as "carbuncles" , etc. 

The states in the cell centers are defined by the conservative difference scheme (|T5|) . 
To obtain the second-order temporal accuracy we must average the wave amplitudes 
over the corresponding domains of influence. To define the fluxes we need to solve the 
Riemann problem between the states at the cell interfaces computed with the averaged 
amplitudes. 

To preserve zero divergence of the magnetic field in three dimensions, we use an 
unsplit Godunov method based on the constrained transport approach [l8[ . We use the 
information about the magnetic field gradients to fulfill the constraint on the magnetic 
field more accurately. 

We tested the performance of PPML on several numerical problems which demon- 
strated its high accuracy on both smooth and discontinuous solutions. Two-dimensional 
flow fields generated by PPML are highly resolved without any wiggled contour lines. Our 
pilot simulations of supersonic magnetized turbulence in three dimensions with PPML 
show that low dissipation and wide spectral bandwidth of this method make it an ideal 
candidate for direct turbulence simulations. 
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